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ABSTRACT 


Internal  filamentary  glass  daircge  caused  by  high  po'  ^r  Q-switched 
pulse  lasers  and  filamentary  trapping  in  liquids  is  analyzed  theoretically 
in  this  report.  Several  models  are  proposed  and  discussed  for  electrostrictively 
driven  acoustic  trapping.  An  analysis  of  Kerr  effect  trapping  is  also  given 
for  purposes  of  comparison. 

In  the  acoustic  theory,  electrostriction  is  the  sound  wave  driving 
force.  The  sound  wave  compressions  cause  focusing  of  the  light  wave  fields. 

The  focused  light  fields  in  turn  cause  stronger  electrostriction  forces.  When 
the  beam  power  is  large  enough  and  the  laser  pulse  duration  is  approximately 
equal  to  the  time  required  for  sound  to  cross  the  unfocused  beam  radius,  the 
trapping  process  runs  until  the  beam  is  focused  to  a  small  radius  limited  by 
diffraction. 


The  theoretical  trapping  thresholds  are  calculated  from  the  laser 
wavelength  and  the  density,  refractive  index,  Young's  modulus,  and  Poisson's 
ratio  of  a  uolid  material,  or  the  density,  refractive  index,  and  speed  of  sound 
of  a  liquid  medium.  These  thresholds  agree  with  experimental  glass  damage 
thresholds  to  within  experimental  error,  and  they  vary  the  same  way  with 
initial  beam  size.  Computer  movies  showing  the  formation  of  strongly  focused 
regions  are  presented.  An  explanation  is  given  for  most  of  the  salient  features 
observed  in  the  damage  phenomenon.  Mathematical  analyses  of  various  features 
of  the  models  are  presented  with  computed  graphs. 
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SECTION  I 

PREFACE 


The  goal  of  the  research  reported  here  is  to  continue  the  develop- 
ment  of  a  theory  of  acoustic  beam  trapping.  The  theory  is  a  possible  explana¬ 
tion  of  the  mechanism  for  internal,  filamentary  glass  damage  by  lasers.  It 
may  also  be  an  explanation  for  some  of  the  beam-trapping  phenomena  observed  in 
liquids  when  they  are  traversed  by  a  high  intensity  laser  pulse. 

We  began  developing  the  theory  in  1965.  Sections  II  and  III  pre¬ 
sent  theoretical  results  obtained  by  August  1966. 

Section  IV  covers  an  analytical  solution  to  the  beam-trapping  equa¬ 
tions  for  the  case  in  which  Kerr  effect  trapping  is  the  dominant  trapping  mech¬ 
anism  and  where  electrostrictively  driven  sound  waves  are  weak  or  absent. 

Section  V  presents  some  of  our  computer  results,  depicting  the 
acoustic  beam-trapping  phenomenon.  An  explanation  of  the  salient  features  of 
the  glass  damage  phenomenon  is  given. 

Section  VI  presents  derivations  of  formulas  used  in  the  computer 
programs.  It  is  the  summary  of  work  performed  l.i  the  second  half  of  the  contract. 
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SECTION  II 

INTRODUCTION 


2.1  A  NEW  POWER  LIMIT 

As  higher  and  higher  power  lasers  are  developed,  basic  limitations 
in  the  power-transmitting  capability  of  materials  and  propagation  media  are  be¬ 
ing  discovered.  Examples  of  limitations  are  electric  breakdow  .  and  beam  insta¬ 
bility.  During  our  study  of  gain  saturation  and  other  anomalies  in  stimulated 
Raman  effect  and  in  our  experimental  work  in  laser  damage  to  glass,  we  have 
identified  a  new  kind  of  beam  instability,  acoustic  beam  trapping. 

Acoustic  beam  trapping  is  caused  by  the  focusing  action  of  electro- 
strictively  driven  acoustic  waves.  These  waves  alter  the  index  of  refraction 
of  the  medium  by  the  density  changes  they  cause.  All  optical  materials  experi¬ 
ence  these  electrostriction  forces  and  photoelastic  effects. 

Acoustic  beam  trapping  sets  upper  limits  to  beam  power  that  can  be 
transmitted  in  important  materials  such  as  glass  and  air.  Furthermore,  the 
scaling  laws  for  acoustic  trapping  (power  threshold  versus  beam  size  and  pulse 
half-time)  differ  from  the  scaling  laws  for  other  types  of  beam  instability, 
such  as  Kerr  effect  or  anomalous  dispersion  trapping.  Ir.  fact,  for  many  trans¬ 
parent  materials,  there  is  a  beam  size  and  pulse  half-time  for  which  the  acoustic 
trapping  threshold  is  lower  than  the  thresholds  for  other  known  instabilities. 

2.2  THRESHOLD  PREDICTABLE 

At  the  present  time,  we  can  predict  the  acoustic  trapping  power 
threshold  in  glassy  materials  from  a  knowledge  of  material  properties,  wave¬ 
length,  pulse  half-time,  and  beam  size.  There  is  a  minimum  power  level 
for  each  material  which  can  be  calculated  from  tabulated  material  properties 
and  the  laser  wavelength. 
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2.3  RECONCILABLE  WITH  STEAD!  STATE  TRAPPING  THEORY 

We  have  also  shown  that  for  any  given  beam  size  and  type  of  material, 
there  5s  a  maximum  average  rate  of  increase  of  power  that  can  be  transmitted 
wicnout  causing  acoustic  beam  trapping.  So  long  as  power  is  a  Jed  to  a  beam  at 
a  slower  rate  than  this,  it  is  oossible  in  principle  to  reach  the  steady  state 
trapping  threshold  predicted  by  Chiao,  Garmire,  and  Townes  (Ref.  1). 

2.4  SCALING  LAWS 

Figure  1  illustrates  tha  scaling  law3  for  three  different  effects 
that  limit  the  ability  of  an  optical  material  to  transmit  a  laser  beam.  In 
the  graph,  total  beam  power  is  plotted  versus  beam  radius  on  log- log  scales. 

The  three  effects  are  electric  field  breakdown,  Kerr  effect  trapping  and  acous¬ 
tic  trapping. 


For  electric  field  breakdown,  as  in  laser- induced  "air  sparks," 

2 

Power  Threshold  _  f  L  nstant  for  *\  f  Beam  'N 
for  Breakdown  -  V.  Material  J  VRadius/ 

On  the  log- log  plot,  this  power  threshold  is  a  line  with  a  slope  of  2.  Thus, 
electric  field  b re V* down  limits  the  material  to  transmission  of  power  levels 
and  beam  sizes  in  the  right-hand  portion  of  the  graph. 


For  Kerr  effect  trapping,  as  shown  by  Chiao,  Garmire,  and  Townes 

(Ref.  1) 


Power  Threshold  for 
Kerr  Effect  Trapping 


(W s  clength) 


(Constant  for  'N 
Material  J 


Since  the  threshold  is  independent  of  beam  size,  it  appears  as  a  horizontal  line 
on  tha  graph.  For  transmission  without  Karr  effect  trapping,  the  beam  power 
level  must  be  below  the  Kerr  effect  threshold. 


Together,  Kerr  effect  trapping  and  electric  field  breakdown  limit 
beam  transmission  to  beam  sizes  and  power  levels  in  the  lower  right  portion  of 
the  graph.  Br  h  of  the  effects  are  virtually  independent  of  the  pulse  half¬ 
time. 
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The  acoustic  trapping  threshold  in  be  calculated  as  in  Section  III 
as  long  as  the  laser  pulse  half-time,  p,  is  less  than  the  time  required  for 
sound  to  cross  the  beam  radius.  The  scaling  law  for  acoustic  trapping  is 


Power  Threshold  for 
Acoustic  Trapping 


2 

(Wavelength)  ( 


Constant  for 
Material 


(Beam  Radius) 
(Pulse  Half-time)^ 


for  the  limiting  case  of  a  ve‘ y  short  pulse  or  a  very  large  beam.  If  the  pulse 
half-time  is  held  constant,  the  threshold  power  can  be  plotted  versus  beam  size 
as  shown  in  the  graph.  The  acoustic  trapping  power  threshold  curve  flattens  out 
at  the  bottom,  at  the  minimum  power  level.  This  minimum  power  level,  P  ,  is 
characteri  only  of  the  material  and  the  wavelength. 


For  the  domain  in  which  the  pulse  half-time  is  longer  than  the  time 
T  required  for  sound  to  cross  the  beam,  the  acoustic  trapping  threshold  remains 
at  P  .  This  power  level  can  be  exceeded  only  if  the  power  does  not  increase 
more  rap illy  than  Fmin/T  in  any  time  T. 

If  the  pulse  half-time  is  increased  by  a  numerical  factor,  and  the 
beam  radius  is  increased  by  the  same  factor,  the  power  threshold  remains  con¬ 
stant.  Thus,  in  the  log-log  graph,  the  curve  is  shifted  to  the  right  by  the 
log  of  the  numerical  factor  by  which  the  beam  radius  is  increased.  Thus  the 
graph  shows  three  extra  curves  shifted  by  factors  of  10,  100,  and  1000. 

In  some  materials  the  Kerr  effect  trapping  threshold  is  lower  than 
the  minimum  power  level  for  acoustic  trapping.  In  these  materials  Kerr  effect 
trapping  will  occur  before  acoustic  trapping.  However,  most  common  optical 
materials  have  a  critical  power  level  lower  than  the  Kerr  effect  trapping  thresh¬ 
old.  For  these  materials  there  will  always  be  a  domain  of  beam  size  and  pulse 
half-time  in  which  acoustic  bean,  trapping  sets  the  maximum  transmittable 
power  for  the  material. 

2.5  TYPICAL  RESULTS  FOR  GLASS 

Acoustic  trapping  is  an  important  cause  of  laser  damage  to  optical 
glass,  as  shown  in  experiments  performed  by  Steinberg,  Atwood,  Lee,  and  Ward 
(Ref.  2).  The  theoretical  trapping  threshold  is  compared  with  the  experimental 


6 


PERKIN-ELMER 


Report  No.  9204 


damage  threshold  in  Figure  2.  Generally,  the  trapping  threshold  is  below  the 
experimental  damage  threshold  points,  as  expected.  In  particular,  note  the 
experimental  damage  threshold  curve  for  dense  flint.  The  left,  lower  portion 
has  a  slope  of  2.  In  that  region,  the  cause  of  damage  is  probab'iy  electric 
breakdown.  The  right,  upper  portion  fits  the  curve  for  acoustic  trapping  to 
within  the  experimental  repeatability;  and  it  scales  the  same  way.  For  fused 
silica  and  BK-7  the  agreement  between  the  acoustic  trapping  threshold  and  the 
measured  damage  threshold  is  even  better. 

The  experimental  results  are  not  attributable  to  Kerr  effect  trap¬ 
ping  because  of  the  dependence  on  beat  size.  Also,  in  glass,  Kerr  effect  is 
so  weak  th3t  the  power  threshold  for  Kerr  effect  trapping  is  4  megawatts,  which 
is  above  the  top  of  the  graph. 
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SECTION  III 


DERIVATION  OF  THE  ACOUSTIC 
BEAM  TRAPPING  THRESHOLD 


3.1  THE  DRIVING  FORCE 

Electrostriction  is  the  force  exerted  by  an  electric  field  on  a 
material  medium,  when  the  force  is  proportional  to  the  square  of  the  field. 
The  net  body  forcu  is  then  proportional  to  the  power  intensity  gradient,  in 
a  lossless  medium.  The  relevant  permittivity  ratio  at  optical  frequencies  is 
the  square  of  the  refractive  index.  Thus,  the  net  body  force  f  per  unit 
volume  due  to  the  light  beam  is  (Ref.  3) 

f  =  (l/6n  c)  (n2  +  2)  (n2  -  1)  “  I. 

o  o  o 

Thus  a  cylindrical  beam  will  d'-ive  a  radially  propagating  sound  wave. 

3.2  THE  PHOTOELASTIC  EFFECT 

The  refractive  index  change  due  to  small  acoustic  compressions 
is  defined  as 


in  =  n  -  n  =  (p  -  p  )  Sn/So 

Assuming  constant  polarizability  per  molecule,  An /Ap  may  be  calculated  by 
differentiating  the  Clausius-Mosotti  relation.  Then  we  obtain 

^n  =  (1/6  nQ)(n2  2)  (n2  -  1)  a 

where  (y  is  the  normalized  compression. 

a  =  (p  -  p0)  /p0 
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Other  effects,  such  as  the  Kerr  effect,  may  add  to  An.  They  will 
not  have  the  same  distribution  as  acoustic  compression,  in  general,  nor  will 
they  vary  the  same  way  with  beam  size.  However,  they  can  be  added  to  An  later. 
Let  us  ignore  them  for  the  present. 

3.3  THE  LASER  BEAM 

For  the  purpose  of  this  discussion,  the  laser  pulse  is  represented 
by  two  impulses  of  equal  energy.  Each  impulse  has  a  Gaussian  radial  intensity 
distribution.  The  radius  is  measured  to  the  point  where  the  intensity  drops 
to  1/e  of  the  peak  intensity.  Tne  beam  is  circularly  symmetrical  and  gently 
focused.  At  the  focus  the  intensity  distribution  is 

I(r,  t)  =  (W/it  w1 2 3 4 5)  exp  (-r2/w2)^6  (t)  +  6(t-At)l  /2 

where  W  is  the  energy  of  two  impulses.  The  two  impulses  are  separated  by  a 
time,  At,  equal  to  2/3  the  half-height  duration  or  half-time,  p,  of  the  laser 
pulse.  Thus  if  the  physical  laser  pulse  has  a  Gaussian  time  distribution,  the 
two  impulses  will  occur  at  the  centroids  of  the  two  halves  of  the  physical  pulse. 

The  temporal  intensity  distribution  is  shown  in  Figure  3,  and  the 
spatial  intensity  distribution  is  sketched  in  Figure  4.  The  resulting  force 
distribution  is  given  in  Figure  5. 

3.4  THE  ACOUSTIC  WAVE 

The  electrostrictive  force  drives  a  radially  propagating  acoustic 
wave.  Usually  when  acoustic  beam  trapping  occurs  the  boundaries  of  the  mate¬ 
rial  are  so  far  from  the  beam  center  that  the  trapping  event  occurs  sooner 
than  sound  can  be  reflected  from  the  boundaries.  Hence,  the  acoustic  wave 
equation  applies  with  the  following  conditions: 

1)  Circular  symmetry 

2)  Infinite  homogeneous  isotropic  medium 

3)  Solution  at  the  beam  center  is  well  behaved 

4)  Solution  at  infinite  radius  is  zero 

5)  Solution  is  not  a  function  of  z  (paraxial  beam  case) 
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Figure  3.  Laser  Pulse  Temporal  Intensity  Distribution, 
and  the  Two-Impulse  Approximation. 
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Figure  4. 

Gaussian  Spatial  Intensity 
Distribution 


Figure  5. 

Initial  Force 
Distribution 


Figure  6. 

Initial  Rate  of  Change 
of  Compression  Distribution 


Figure  7, 

Hankel  Transform  of 
Initial  Rate  of  Change 
of  Compression 


Figure  8. 

Acoustic  Compression  at 
T  =  0.596.  Proportional 
to  Refractive  Index  Change 
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The  acoustic  wave  equation  can  be  solved  for  the  case  of  an  impulse 
driving  force.  Since  the  equation  is  linear,  it  is  convenient  to  use  dimen¬ 
sionless  variables  and  a  normalized  impulse. 

Let  the  dimensionless  radial  coordinate  be  x,  the  radius  measured 
in  units  of  the  initial  focused  beam  radius  w; 

x  =  r  /w 

There  is  a  characteristic  acoustic  response  time.  T,  equal  to  the  time  re¬ 
quired  for  sound  to  traverse  the  beam  radius; 

T  =  w/v 

where  v  is  the  speed  of  sound  for  a  two-dimensional  compression  wave. 

v  =  (Y/2 p)in  (1  -  c)'1/4  (1) 

Here  Y  is  Young's  modulus  and  g  is  Poisson's  ratio  for  solid  media. 

'“t.  the  dimensionless  time  variable  be  r,  time  measured  in  units  of  T. 

T  =  t/T  =  (v/w)  t 

This  choice  of  units  makes  the  velocity  of  sound  equal  unity  in 
the  wave  equation.  Also,  the  acoustic  velocity  u  is  normalized  by  the  velocity 
of  sound,  and  the  acoustic  displacement  u  is  normalized  by  w.  The  acoustic 
compression  CT  is  the  negative  divergence  of  the  displacement: 

O  -  -  V  •  TJ 

where  U  =  u/w.  Differentiating  with  respect  to  r,  we  have 
3a/dT  -  -  V  •  SU /r)T 

where  oU/^r  =  u/v 


13 


PERKirJ-ELMER 


« 


Report  No.  9204 


The  initial  rate  of  change  of  compression  distribution  appears  in 

Figure  6. 


When  the  electrostrictive  force  acts  impulsively  on  a  medium  ini¬ 
tially  at  rest,  there  is  no  immediate  displacement  or  compression.  The  ini¬ 
tial  velocity  and  initial  rate  of  change  of  compression  may  be  deduced  from 

conservation  of  momentum, 
o+ 


'  f  dt 


(x,0+) 


The  compression  a  is  a  function  of  radius  and  time  only.  To  sum 
up,  the  problem  is  to  solve  the  acoustic  wave  equation; 


2 

v  a 


2  , 

d  a/9v 


2 


subject  to 

1)  a  =  a(x,  t) 

2)  a(x,  o)  =  0 

3)  da(x,0)/3r  =  A(l  ■  x2)  exP  (-x2) 

where 

A  =  (  W/3nn  cvp  w^)  (n2  +  2)  (n^  -  1) 

o  o  o  o 

4)  a  (o,t)  4  o c 

5)  ate ,  t)  =  o 

The  equation  and  all  of  the  conditions  except  (3)  are  satisfied 
by  a  (x,  t)  =  (B/A)J  (yx)  sin  (yr)  where  B  and  y  are  arbitrary,,  A  linear 
superposition  of  solutions  will  satisfy  condition  (3)  as  well.  The  correct 
combination  of  solutions  g(y)  is  given  by  the  Hankel  transform  (not  the 
Fourier  transform  because  of  circular  symmetry).  The  weight  function  is  x 
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and  the  eigenfunctions  are  J  (yx)  corresponding  to  1  and  sin  (xy)  or  cos  (vx) 

o 

for  a  Fourier  transform.  Thus, 

or 

8<y)  =  F  <■  ["A(l  -  E?)  exp  (-£2)  1  J  <£y)d£ 

O 

2  2 
-  A(y  /6)  exp  (-y  /4) 

This  distribution  is  plotted  in  Figure  7.  Now  superposition  is 
applied  to  obtain  the  solution. 

(*  **  2  2 

CT(x,r)  =  A  I  (y  /8)  exp  (-y  /4)  J  (yx)  sin  (y r)  dy 
o 

This  difficult  integral  has  been  evaluated  numerically  as  described  below. 

For  the  important  on-axis  case,  where  x  =  0,  the  integration  may  be  performed 
analytically,  since  JQ(^)  =  1* 

80 

(7(0,  t)  =  A  f  (y2/8)  exp  (-y2/4)  sin  (yr)  dy 
% 

=  A  r/2  +  (A/2)  (1  -  2r2)  exp  (-T2)  !*  ^xp  (£2  )  d£ 

’  o 

=  A  (r-  4r'/3  +-•--) 

This  function  is  plotted  in  Figure  9. 

3.5  THE  NUMERICAL  INTEGRATION 

The  integral  for  tf(x,  y)  was  computed  in  our  Scientific  Computer 
Facility.  We  programmed  the  Scientific  Data  Systems  9300  Computer  with  65 

Fortran  IV  statements.  Simpson's  method  was  used.  The  Bessel  function  was 

_8 

generated  by  a  polynomial  approximation  accurate  to  5  x  10  absolute  error 
from  the  Handbook  of  Mathematical  Functions  (Ref.  4,rp.j  3C9f).  One  hundred 
and  ona  values  of  y  were  used  each  time  the  integral  was  evalutated.  The  in¬ 
tegral  was  computed  et  101  values  of  x  for  each  of  124  values  of  t.  Time  was 
3aved  by  storing  parts  of  the  kernel  that  do  not  change,  and  by  buffering  the 
printer.  Running  time  was  about  40> minutes,  or  about  1  millisecond  for  each 
evaluation  of  the  Bessel  function. 


15 


PERKiN-ELMER 

Report  No.  9204 


The  computed  points  were  saved  on  magnetic  tape.  We  displayed  the 

points  wi  h  the  on-line  oscilloscope  associated  with  our  SDS  930  Computer. 

The  Fortran  II  display  program  allowed  us  to  show  stationary  frames  continu¬ 
ously,  or  to  show  frames  sequentially  at  any  comfortable  viewing  rate.  The 
display  was  useful  in  studying  the  wave  motion  and  in  checking  the  computati  n. 

Refer  to  Figures  10  and  11  for  plots  of  wave  amplitude  a/A  versus 
x  at  various  values  of  r.  The  compression  at  maximum  on-axis  amplitude  is 
plotted  in  Figure  8.  for  comparison  with  the  other  functions  in  the  solution. 

3.6  THE  TRAPPING  CONDITION 

An  unexpected  result  of  the  computed  solution  of  the  acoustic  wave 
is  the  fact  that  the  compression  is  greatest  on  axis,  for  0  r  ^  0.85.  Thus, 

it  is,  from  the  start,  a  focusing  distribution.  There  is  no  latency  period  of 

zero  or  negative  focusing  before  some  positive  focusing  begins.  The  on-axi.; 
solution  shown  in  Figure  9  shows  that  the  convergence  (reciprocal  focal  length) 
varies  linearly  with  j,  for  j  «  1. 

The  two-impulse  model  is  valid  for  T  <  0.85.  The  first  impulse 
starts  an  acoustic  wave,  and  the  acoustic  wave  builds  up  a  distributed  lens 
<-,f  increasing  convergence.  Diffraction  of  the  second  impulse  will  be  defeated 
if,  at  the  time  the  second  impulse  occurs, 

in/n  >  (1/2)  (\/2^nQw)  2. 

This  on-axis  in  will  be  the  sum  of  in  caused  by  the  acoustic  wave  and  Kerr 
effect.  The  acoustic  in  is  a  function  of  beam  energy,  pulse  duration,  and  ini¬ 
tial  focused  spot  size.  Kerr  effect  depends  on  power.  Thus  the  two  effects  do 
not  vary  the  same  way  with  beam  size. 

3.7  THE  TRAPPING  THRESHOLD 

For  optical  glass  and  o:her  materials  in  which  Keir  effect  is  weak, 
we  may  calculate  the  approximate  energy  threshold  for  trapping  by  equating  the 
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Figure  10.  Sound  Wave  Compression  at  Various  Times 
After  an  Initial  Impulse 
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Figure  11.  Sound  Wave  Compression  at  Various  Times 
Near  the  First  Peak 
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required  axial  an  with  i_he  An  due  to  the  acoustic  wave. 

„  _  9c\2  r  _ »o  _  "|  yw 

TRAP  L  .  2  .2,  2  ..2  o/A 

(nQ  +  2)  (nQ  -  1.) 

We  may  define  the  power  for  the  impulse  pair  as 

P  «  W/at  =  3W/2p  - 

and  substitute  T  *  w/v  so  the  trapping  power  becomes 

Q  .  2  _  n  P  2 

p  _  j  _ °  o _  ]  v _ 

TRAP  4ji  L  ,  2  ,.2,  2  ..2  1  (2P/3T)(o/a) 

(no  +  2)  (nQ  -  1) 

The  minimum  trapping  power  will  be  required  when  the  pulse  half-time  is  such 
that  (2p/3T) (o/A^  is  a  maximum.  This  maximum  is  .261,  and  it  occurs  when  p  =  1.22T 

P  tn  =  8.64c\2  n  p  v2/\n(.n  +  2)2(n2  -  l)2  1 
MIN  oo  L  o  o  l 

This  minimum  power  depends  only  on  the  material  properties  and  the  wavelength. 

The  beam  size  must  be  matched  to  the  pulse  half-time  as  follows: 

p  =  1.22T  =  1.22  w/v 

w  ,  .  =  .82  pv 

matched  r 

For  values  of  p/T  «1,  o/A  —  2p/3T. 

Hence  in  the  limit  of  short  pulses  or  large  beams, 

t,  2.2 

PTRAP  varies  as  w  /?  • 
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3.8  NUMERICAL  CALCULATIONS 

The  acoustic  trapping  threshold  curve  for  a  material  ran  be  ob¬ 
tained  with  simple  calculations.  It  is  only  necessary  to  obtain  the  sound 
velocity,  the  beam  size  that  matches  the  pulse  half-time,  and  the  minimum 
power  level.  Then  simple  graphical  methods  yield  the  threshold  curve. 


Sound  Velocity:  The  relevant  sound  velocity  is  that  for  a  radially 
propagating  compressional  wave.  For  solids,  use  equation  (1)  in  Section  3.4. 
The  Young's  modulus,  density,  and  Poisson's  ratio  ror  glasses  are  found  in  the 
newer  Schott  glass  catalogs. 

2  3 

Example:  For  BK-7,  Y  =  8310  kp/mm  ,  p  =  2.51  g/cm  ,  f  =  0.208. 

3 

Multiply  p  by  10  to  obtain  the  density  in  kilograms  per  cubic  meter,  and  Y 

by  9.81  x  10^  to  obtain  newtons  per  square  meter  from  kiloponds  per  square 

3 

millimeter.  Hence  v  =  4.27  x  10  meters  per  second. 


Matched  Beam  Size:  Let  the  pulse  half-time  p  be  measured  at 


"matched  Pv' 


the  half-peak  points.  Then  w^ 

Minimum  Power  Level  The  minimum  power  level  for  the  material 

2  1 


is 


PMIN  =  <8-6Acx2/*>  Pvn0/[(no+2)  (nc 


1) 


Example:  For  BK-7,  at  \  =  694.3  nm,  n^  =  1.45.  The  critical 
power  level  is  365  kilowatts. 


Plotting  the  Threshold:  Use  log- log  graph  paper  having  equal  size 
divisions  for  the  horizontal  and  vertical  scales.  Draw  a  horizontal  line 
corresponding  to  the  minimum  power  level.  Mark  the  matched  point  on  the 
minimum  power  line,  corresponding  to  the  matched  beam  size. 


For  a  quick,  conservative  overestimate  of  the  trapping  threshold 
curve,  simply  draw  a  line  with  a  slope  of  2  upward  from  the  matched  point. 
(The  angle  is  63.5°  to  the  horizontal.) 
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For  an  accurate  graph,  plot  the  curve  in  Table  I  on  a  second 
piece  of  log-log  paper.  Lay  the  first  piece  over  the  second  so  the  axes  are 
parallel  and  the  matched  point  is  just  over  che  bottom  point  on  the  second 
graph.  Trace  the  curve. 


TABLE  I 

NUMERICAL  VALUES  USEFUL  FOR 
PLOTTING  ACOUSTIC  TRAPPING  THRESHOLD  CURVES 


Horizontal  Axis 

0.91 

1.00 

1.10 

1.25 

1.43 

1.68 

2.00 

2.50 

3.33 

5.00 

10.0 


Vertical  Axis 

0.965 
1.00 
1.05 
1.20 
1.43 
1.78 
"  .  36 
3.48 
5.90 
12.8 
55.0 
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SECTION  IV 

AN  ANALYTICAL  SOLUTION  OF  THE  LASER  BEAK 
RADIUS  EQUATION  FOR  KERR  EFFECT  TRAPPING 


A  high  intensity  laser  beam  passing  through  a  material  medium  can 
focus  itself  into  a  long,  thin  filament  and  propagate  without  normal  diffrac¬ 
tion  spreading.  This  self-trapping  phenomenon  arises  when  the  medium's  re¬ 
fractive  index  n  is  higher  along  the  beam  axis  than  along  the  beam  edges. 

Such  refractive  index  distribution  acts  like  a  series  of  thin  positive  lenses, 
o,-?  shown  in  Figure  12. 

The  high  intensity  beam  can  set  up  a  focusing  refractive  index 
distribution  by  several  physical  mechanisms,  such  as  electrostriction,  anom¬ 
alous  dispersion,  and  reorientation  of  molecular  dipol?  moments..  The  last 
effect  is  Kerr  effect,  after  its  discoverer,  John  Kerr  (1824-1907).  In 
liquids  the  molecular  reorientation  can  occur  in  times  on  the  order  of  10 
picoseconds.  The  effect  is  thus  virtually  instantaneous  compared  with  the 
duration  of  nanosecond  laser  pulses,  but  it  is  much  too  slow  to  follow  500 
terahertz  light  wave  fields.  The  local  change  in  index  of  refraction  is 
proportional  to  the  local  beam  intensity.  When  the  beam  power  is  above  a 
certain  threshold  power  level,  the  laser  beam  focuses  itself  to  a  smaller 
beam  size  and  higher  intensity.  The  smaller,  higher  intensity  beam  causes 
still  stronger  self-focusing  until  the  beam  has  trapped  itself  at  a  small 
radius,  limited  by  diffraction. 

4.1  THE  BEAM  RADIUS  EQUATION 

This  note  sets  forth  an  analysis  of  Kerr  effect  trapping,  The 
analysis  is  based  on  a  ray  tracing  equation  published  by  Tien,  Gordon,  and 
Whinnery  (Ref.  6).  Equation  (9)  in  that  article  reduces  to 

2  2  -3  2  2  > 

9  a/^z  =  a  (\/2-rno)  +  (a/n^  n(0,  z,  t) /?*r“  (2) 


23 


Equivalent  Beam  Guiles.  The  equivalence  may  be  seen  by  integrating  the  optical 
path  length  (n  dz )  along  the  axes  and  along  the  edges  of  both  figures.  In  both 
there  is  more  retardation  on  axis  than  off. 
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where 


a(z,t) 

n(r,z,t) 

n 

o 

v. 

r,z 


beam  radius  to  the  point  where  the  intensity  is 
1/e  of  the  peak  intensity, 
local  index  of  refraction 

undisturbed  or  nominal  index  of  the  medium 
vacuum  wavelength  of  the  laser  beam, 
cylindrical  coordinates  for  the  beam. 


For  this  equation  to  be  valid  Tien,  Gordon,  and  Whinnery  require 

that 


1)  "the  light  beam  is. .. launched  with  a  Laguerre-Gaussian 
or  Hermife-Gaussian  field  distribution" 

2)  "the  refractive  index  of  the  medium  varies  slowly  in  space 
(negligibly  in  an  optical  wavelength)" 

3)  the  variations  in  n  are  small  compared  with  n 

o 

The  last  two  restrictions  are  violated  when  the  beam  is  trapped 
to  a  filament  with  approximately  one  wavelength  radius,  and  the  scattered 
light  is  spread  over  a  wide  wavelength  range.  However,  the  equation  can  be 
used  to  study  the  collapse  of  the  beam  toward  the  trapped  condition,  and  to 
study  the  phenomenon  when  the  beam  power  is  below  threshold. 


4.2  BEAM  PROPAGATION  WITHOUT  SELF-FOCUSING 


When  the  beam  intensity  is  weak,  the  last  term  in  equation  (2) 
is  negligible.  The  ray  paths  are  then  hyperbolas.  The  beam  radius  solution 
is 


1/2 


a(z)  =  |^(a  + 


a1“) 


(\z/2itn  a  )' 
o  o 


(3) 


where  a  and  a.  are,  respectively,  the 
o  1 

rays  just  after  entering  the  medium  at 

n  =  1.5,  a  =  0.1mm,  a  =  -0.001839. 
o  o  1 

0.05mm  at  z  -=  4.1cm.  This  is  the  case 


initial  radius  and  initial  slope  of  the 
z  =  0 .  For  example,  let  \  =  lum, 

Then  the  beam  is  focused  to  a  radius  of 
of  a  very  gently  focused  beam.  The 
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first  frame  of  Figure  14  shows  a  plot  of  a(z)  for  this  case  .  Note  that  the 
vertical,  radius  scale  is  greatly  exaggerated. 

4.3  THE  LASER  BEAM  INTENSITY  EQUATION 

The  lowest  order  mode  for  the  laser  beam  has  a  nearly  Gaussian 
intensity  profile  characterized  by  the  radius,  a,  to  the  point  where  the 
intensity  drops  to  1/e  of  the  peak  intensity.  If  the  propagation  medium  is 
not  too  violently  inhomogeneous  and  if  the  inhomogeneity  is  radially  sym¬ 
metrical,  the  Gaussian  profile  is  maintained  along  the  beam,  although  the 
beam  radius  varies  because  of  diffraction  and  focusing  by  the  inhomogeneities. 

The  laser  power  may  also  vary  as  a  function  of  time.  For  Q-switched 
pulse  lasers  generally  the  pulse  energy  can  be  measured  and  some  rough  idea  of 
the  time  distribution  of  the  pulse  can  be  obtained.  For  purposes  of  analysis 
we  can  use  a  simple  pulse  shape  such  as  the  unit  quartic  pulse  shown  in 
Figure  13.  This  shape  has  a  continuous  derivative  and  finite  extent. 


Thus  the  laser  beam  intensity  is 

I(r, z, t)  =  (W/na2)  exp  (-r2/a2)  q  ( t ; p) 


where  W  is  the  pulse  energy  and  q ( t ; p)  is  the  unit  quartic  pulse  of  time 
constant  p: 


q(t;p) 


C 


(15/16p)  (t/p) 2 ( 2- t / p) 2,  0  £  t  s  2p 
0  otherwise 


Of  course,  a  is  a(z,t)  given  by  the  solution  of  equation  (2).  The  intensity 
distribution  is  so  normalized  that 


r 2rt 

d9 


r  r  dr 

-o 


dt  I  (r , z, t)  =  W 


independent  of  the  value  of  z. 
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i/p  q(t;p) 


/  (15/16p)  (t/p)2(2-t/p)2  for  0  <r  t  <  2p 
1  0  otherwise 

f  9(t;p)  dt  =  1 

J  -to 


Figure  13.  Unit  Ouartic  Pulse 
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4,4  VARIATION  OF  REFRACTIVE  INDEX  WITH  INTENSITY 

The  change  in  refractive  index  can  be  expressed  in  cgs  units  as 
2  ,„2 

n  =  n  +  ~  \JE  +  ..... 

where  J  is  the  high  frequency  Kerr  constant  due  to  molecular  rotation.  Thus 
the  change  in  index  is  proportional  to  the  local  intensity  for  either  a  plane 
polarized  or  circularly  polarized  oeam  (although  the  constants  of  proportion¬ 
ality  differ  for  the  two  case?).  Let  the  constant  of  proportionality  be  K, 
depending  on  the  wavelength,  the  material,  and  the  beam  polarization. 

n(r,  z,t)  =  nQ  +  KI  (r,z,t) 

In  equation  (2)  the  last  term  Becomes 
-a"3(2  KW/rtnQ)  q  (t;p) 

Since  both  tei ms  on  the  right  of  equation  (2)  are  proportional  to 
-3 

a  ,  the  solution  at  any  time  will  be  a  hyperbola  whose  value  and  slope  at 
z-.  0  are  a^  and  a^,  respectively.  However,  the  focal  point  and  semiaxes  of  the 
hyperbola  may  vary  with  time.  The  equation  becomes 

2  2  -3 

d  a/*z  =  f (t) a 

where 

f(t)  =  (\/2nno)2  -  2KWq(t;p)/nno. 

The  solution  is 

1/2 

T  2  -2  2  1  ' 

a (z, t)  =  I  (a  +  a  z)  +  f(t)  a  z 

C(  I  o 

A  negative,  zero  or  complex  value  of  a  would  be  a  r.on-physical 

solution.  This  places  a  restriction  on  the  minimum  value  of  f(t)  which  occurs 

2 

when  t  =  p.  We  have  a  (z,p)  equal  to  a  quadratic  expression  i.  z.  In  order 

that  a^(z,p)  be  greater  than  zr.ro  for  all  real  values  of  u,  the  discriminant 

2 

of  the  quadratic,  expressior  must  be  negative.  Then  the  only  zeros  of  a  will 
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be  complex.  This  condition  implies  th^t  f(t)  must  always  be  positive.  The 
minimum  value  of  f(t)  occurs  when  t=p,  i.e.,  at  the  peak  of  the  laser  pulse. 
For  a  valid  solution  we  must  have 


/n/n  <  -y  (\/2im  a)  2 
o  l  o 

where  ^n,  the  on-axis  increase  in  refractive  index,  is 
c.n  =  KWq(p;p)  /.t a"\ 


This  condition  is  equivalent  to  the  Pierce  stability  criterion  L/f  >'  4.  where 
L  and  f  are  the  spacing  and  focal  length  of  lenses  in  a  sequence  (Ref.  5). 

2 

The  threshold  value  of  KW/p  becomes  (KW/p)  =  2\  /tn 

threshold 


Let  k  be  the  fraction 


k  = 


(KW/p) /(KW/p) 


threshold 


Then  f(t)  becomes 

f(t)  ---  (\/2:no)2  [l-kq(t;p)/q(p;p)3 

Refer  to  Figure.'  14  through  16  for  plots  of  the  beam  radius  when  k  =  0.999. 
0.95,  an’  0.85. 
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SECTION  V 

BEAM  TRACING  IN  THE  ACOUSTIC  TRAP 

The  theory,  developed  from  the  two-impulse  model  in  Sections  It 
and  III.  yields  a  threshold  for  acoustic  beam  trapping  which  is  below  the  ex¬ 
perimental  damage  thresholds  for  the  glassy  materials  tested.  It  was  solved 
by  completely  analytical  means,  using  tabulated  properties  of  special  functions. 
The  computer  was  only  used  to  obtain  values  for  the  radial  sound  wave  off  the 
axis  of  the  beam.  That  computation  was  not  essential  to  the  solution,  although 
it  did  give  the  following  important  insight:  The  refractive  index  distribution 
set  up  by  the  acoustic  wave  is  initially  a  focusing  distribution,  although  its 
strength  varies  with  time.  However,  the  theory  did  nor  allow  for  the  continuous 
interaction  of  light  and  sound,  nor  did  it  show  how  the  focusing  action  develops 
with  time. 

5.1  OTHER  FEATURES  OF  THE  GLASS  DAMAGE  PHENOMENON 

A  complete  theory  shou’d  also  explain  the  other  observed  features 
of  internal,  filamentary  glass  damage,  such  as  the  following: 

1)  The  Spectrum  of  the  Side-Scattered  Light 

A  white  light  flash  is  seen  when  the  damage  event  occurs. 

u  o 

Ruby  laser  light  may  be  shifted  from  6943A  down  to  4000A. 

Such  large  shifts  can  only  be  explained  by  a  very  strong 
dynamic  refractive  index  change. 

2)  Starting  Location  of  the  Damage  Track 

Even  if  the  laser  is  focused  at  the  entrance  face  of  the 
gla  s  sample,  the  damage  track  never  begins  at  the  entrance 
face.  There  is  always  a  short  interval  between  the  entrance 
face  of  the  sample  and  the  start  of  the  track. 
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3)  Exit  Face  Pitting 

Usually  the  track  ends  in  a  pit  on  the  exit  face  of  the 
glass  sample.  As  a  matter  of  fact,  the  threshold  for  the 
exit  face  pitting  appears  to  be  a  little  lower  than  the 
threshold  for  track  formation.  Nevertheless,  the  phenomena 
is  different  from  the  type  of  surface  damage  reported  by 
others. 

4)  Location  of  Damage  Stars 

The  damage  track  is  not  always  continuous;  it  may  start  and 
stop  several  times  in  the  sample.  Often  there  is  a  damage 
star  on  the  upstream  end  of  the  track.  These  damage  stars 
are  localized  regions  of  gross  fracture.  Occasionally 
they  show  discoloration,  indicating  possible  chemical  de¬ 
composition.  If  the  damage  track  extends  upstream  of  the 
damage  star  it  usually  only  does  so  for  a  very  short  distance 
compared  to  the  extension  downstream  from  the  damage  star. 

The  location  of  the  damage  star  is  intriguing.  If  the 
damage  starts  in  the  damage  star  and  then  propagates  down¬ 
stream  to  form  the  track,  why  doesn't  the  damage  star  cast 
a  downstream  shadow  and  prevent  the  beam  from  concentrating 
in  the  thin  filament? 

5)  Track  Propagation  Speed 

Since  the  damage  event  occurs  in  nanoseconds,  its  dynamics 
are  difficult  to  follow.  There  is  experimental  evidence 
that  the  event  that  forms  the  track  propagates  at  about 
10  times  the  speed  of  sound  in  the  glass.  The  experimental 
evidence  does  not  show  unequivocally  that  the  propagation 
is  either  upstream  or  downstream. 

Most  of  these  features  are  explained  qualitatively  in  the  theory  of 
electrostrictively  driven  acoustic  beam  trapping. 
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5.2  DEFINITION  OF  THE  PROBLEM 

Our  purpose  is  to  understand  how  the  acoustic  beam  trap  forms,  not 
how  the  damage  occurs.  Therefore  we  only  need  consider  small  acoustic  pressure 
changes  in  the  material.  We  will  see  how  these  lead  to  a  beam  instability  when 
the  power  exceeds  a  certain  threshold.  The  nature  of  the  instability  is  such 
that  the  beam  rapidly  focuses  itself  into  a  thin  filament.  The  damage  occurs 
when  the  beam  is  in  this  trapped  condition.  The  instability  itself  can  be 
studied  by  means  of  a  geometric  ray  tracing  equation  which  includes  first  order 
diffraction  effects.  This  beam  tracing  equation  is  equation  (2)  of  Section  IV. 

The  sound  wave  satisfies  essentially  the  same  conditions  as  those 
given  in  Section  III,  except  that  it  is  driven  continuously  by  the  light  wave, 
rather  than  running  inertial ly  after  an  impulse.  For  a  gently  focused  beam, 
axial  components  of  the  sound  wave  field  are  negligible  in  comparison  with 
the  radial  components. 

The  desired  solution  will  show  a  graph  of  the  beam  radius  plotted 
versus  axial  length  for  each  of  many  small  time  intervals  during  the  laser 
pulse.  The  time  distribution  of  the  laser  pulse  can  be  assumed  to  be  the 
quartic  pulse  distribution  described  in  Section  IV,  since  the  actual  details 
of  the  shape  of  the  laser  pulse  do  not  have  too  mu'-h  bearing  on  the  threshold 
for  damage  and  on  the  time  development  ^  -_ap. 

5.3  THE  COUPLED  EQUATIONS 

The  sound  wave  equation  is 

2  2  -1  2  2  -2  2  2 
3  cr/^r  +  r  dor/dr  +  o  ( j/32  =  v  3  (j/^t 

-  ( 1  ^vn^)  (n^  +  2)  (n2  -  1)  (W/it) q( t ; p) 

•  4  a  ^(l-r2/a2)  exp  (-r2/a2) 

The  beam  tracing  equation  is 

32a/3z2  =  a  3(\/27mo>2  +  (a/6n2)  (n2  +  2)  (n2  -  1)  ?j2a(o,  z,  t)  /^r2 
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The  general  solution  of  these  coupled,  nonlinear,  partial  dif¬ 
ferential  equations  has  been  left  as  an  exercise  for  the  computer.  One 
simplification  can  be  made  immediately.  The  beam  tracing  equation  involves 
only  the  second  derivative  with  respect  to  r  of  the  on-axis  sound  wave  solution. 
Since  the  sound  wave  equation  is  linear  we  may  take  its  Hankel  transform  ana¬ 
lytically.  When  z  is  held  fixed,  the  Hankel  transform  of  g  satisfies  an 
ordinary  differential  equation  in  t.  This  ordinary  equation  may  be  solved  using 
initial  values  of  the  Hankel  transform,  and  of  its  derivative  with  respect  to 
time,  and  assuming  a  step-plus-ramp  driving  function.  The  solution  is  then 
found  by  taking  the  inverse  Hankel  transform.  Since  the  beam  tracing  equation 
involves  only  the  second  derivative  of  the  sound  wave  amplitude  at  r  =  0,  and 
the  inverse  Hankel  transform  involves  r  only  as  part  of  the  argument  of  the 
zeroth-order  Bessel  function,  we  may  pt'-form  the  differentiation  analytically 
under  the  integral  sign.  Thus  for  each  time  step  and  each  value  of  z,  the 
sound  wave  equation  can  be  solved  with  only  one  numerical  integration.  The 
use  of  a  step-plus-ramp  driving  function  allows  us  to  take  a  much  coarser 
time  step  than  we  would  require  if  we  used  an  impulsively  driven  equation. 

5.4  THE  STEP-BY-STEP  COMPUTER  SOLUTION 

The  problem  is  divided  up  into  sections  and  slices  as  shown  in 
Figure  17.  In  its  present  state  of  development  the  computer  program  can  handle 
101  sections.  The  zeroth  section  is  the  entrance  face  of  the  solid  or  liquid 
medium. .  The  initial  beam  radius  and  entrance  angle  are  fixed  on  the  zeroth 
section.  Since  the  solution  is  radially  symmetrical,  it  will  be  the  same  in 
each  radial  slice. 

The  computer  program  must  store  values  of  a  for  each  of  the  101 
sections.  Also  for  each  section,  it  must  store  40  values  of  the  Hankel  trans¬ 
form  of  the  sound  wave,  and  40  values  of  the  derivative  of  the  Hankel  transform 
of  the  sound  wave  with  respect  to  time.  An  array  of  101  points  is  also  set 
aside  for  the  second  derivative  of  the  sound  wave  with  rrspect  to  r.  No 
Bessel  function  storage  or  Bessel  function  subprogram  is  required.  Initially 
the  Hankel  transform  of  the  sound  wave  and  its  derivative  with  respect  to  time 
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are  set  to  zero  along  with  the  second  derivative  of  the  sound  wave  with  re¬ 
spect  to  r.  The  solution  then  proceeds  in  steps: 

1)  Solve  the  beam  tracing  equation  using  the  Runge-Kutta  method. 
The  axial  step  size  in  the  solution  is  1/5  of  the  interval 
between  sections.  Values  of  the  second  derivative  of  the 
sound  wave  with  respect  to  r  are  interpolated  by  Lagrange 
cubic  fitting. 

2)  Store  101  values  of  the  beam  radius  a.  Also  save  these 
values  on  tape  for  later  plotting. 

3)  If  the  time  has  not  yet  reached  2. Op,  the  end  of  the 
laser  pulse,  advance  the  time  by  a  small  time  step  and 
obtain  the  new  value  of  the  beam  power  at  that  time. 

4)  For  each  section  of  the  beam  compute  the  second 
derivative  of  the  sound  wave  with  respect  to  r  by 
using  Simpson  intergration.  At  the  same  time  revise 
the  value  of  the  Hankel  transform  of  the  sound  wave 
and  its  time  derivative  for  the  next  time  step. 

5)  Go  to  1 . 

The  program  parameters  are  the  wavelength,  the  nominal  refractive 
ir  !ex,  the  initial  beam  radius  and  the  initial  slope  of  the  beam  radius  at 
the  entrance  face,  and  the  fraction  of  pulse  power  divided  by  theoretical 
threshold. 

5.5  TYPICAL  RESULTS  FOR  ACOUSTIC  BEAM  TRAPPING 

Figure  18  depicts  acoustic  beam  trapping  when  the  pulse  power  is 
3007.  of  the  threshold.  There  are  11  frames  showing  the  beam  radius  plotted 
versus  axial  length  for  11  different  times  during  the  pulse.  The  first  frame 
shows  the  path  of  the  beam  when  the  illumination  has  iust  begun.  The  sixth 
frame  gives  the  trace  when  the  pulse  power  is  highest.  The  11th  frame  shows 
the  path  taken  by  the  light  in  the  trailing  end  of  the  pulse. 


38 


PERKIN-EUtfcER 


Report  No.  9204 


The  action  may  be  interpreted  by  considering  that  each  section  of 
the  beam  forms  an  acoustic  lens.  Initially  the  lenses  have  zero  power,  but 
sone  of  them  build  up  focusing  power  faster  than  others.  The  strongest  lens 
is  formed  a*-  the  focus  of  the  beam.  As  its  focusing  power  builds  up  in  time 
it  causes  the  beam  to  be  refocused  at  some  distance  downstream.  As  the 
strength  of  the  acoustic  lens  continues  to  build  up,  its  focal  length  is 
shortened  so  the  second  focus  moves  upstream.  The  second  focus  likewise  forms 
a  third  focus  and  so  on.  A  little  after  the  laser  pulse  is  half  over,  three 
definite  foci  hav.  appeared  in  the  12-l/2cm  region  plotted.  At  the  same  time 
that  these  acoustic  lenses  are  becoming  stronger  because  of  the  inertial 
properties  of  the  material,  the  laser  pulse  power  is  decreasing.  The  net 
result  is  regions  of  sharper  av.H  sharper  focus  interspersed  by  regions  in  which 
the  beam  radius  is  increasing.  Eventually  the  beam  will  be  able  to  escape 
again  to  large  radii. 

Note  in  the  last  frame  that  the  sharpest  focus  occurs  the  farthest 
downstream.  This  is  to  he  expected  since  the  focus  there  is  the  result  of 
two  strong  lenses  upstream.  Another  important  feature  is  the  motion  of  the 
first  locus.  Notice  how  it  "oves  upstream  in  the  beam.  This  is  because  the 
acoustic  lenses  in  the  first  few  sections  of  the  sample  are  also  developing 
in  strength  with  time. 

Sevei-1  features  of  the  damage  phenomenon  are  inroediatf.ly  ex¬ 
plain  d.  The  fact  that  white  light  is  scattered  largely  from  the  region  of 
the  damage  stars  is  easily  seen  to  follow  from  the  fact  that  the  highest 
focusing  one!  therefore  the  highest  index  of  refraction  occur  at  places  that 
become  damage  stars.  Also  the  track  propagates  backward,  becoming  more  and 
more  sharply  defined  until  it  ends  in  a  dosage  star. -  If  there  are  more  than 
one,  the  damage  stars  occur  in  backward  sequence.  Th’  farthest  downstream 
occurs  first,  then  the  next  one  upstream,  then  finally  the  damage  star  closest 
to  the  upstream  end  of  the  beam.  Thus,  the  damage  e''ent  does  not  cast  a  shadow 
that  gets  tr.  its  own  way. 
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The  reason  why  it  is  difficult  to  measure  the  speed  of  propagation 
of  the  damage  event  is  now  apparent.  The  intensity  varies  widely  and  more  than 
one  focus  region  may  move  past  the  velocity-sensing  optics. 

The  pit  on  the  exit  face  of  the  glass  sample  is  formed  by  axial 
sound  wave  components  which  do  become  large  when  the  beam  is  sharply  focused. 
There  can  never  be  a  damage  pit  on  the  entrance  face.  The  entrance  conditions 
of  the  beam  are  fixed;  therefore,  the  bearr  cannot  focus  itself  sharply  at  the 
entrance  face  of  the  glass. 

All  of  the  qualitative  features  of  the  phenomenon  are  a  sharp  con¬ 
trast  to  the  phenomenon  of  Kerr  effect  trapping,  in  which  there  is  one  focus 
region  that  moves  downstream. 
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SECTION  VI 

DERIVATIONS 


A  number  of  derivations  are  presented  in  this  section.  Each  one 
is  the  analytical  solution  of  some  part  of  the  model  for  acoustic  laser  beam 
trapping.  In  many  cases  the  results  were  evaluated  by  computer  programs.  In 
those  cases,  an  outline  of  the  program  and  typical  plots  are  given. 

6.1  RADIAL  ACOUSTIC  WAVES  UNDER  IMPULSE  EXCITATION  AND  UNDER  CONTINUOUS 
DRIVING  FORCES 

The  crude  (but  useful)  model  of  Section  III  divides  the  laser  pulse 
into  two  impulses.  The  first  impulse  excites  an  acoustic  wave.  By  the  time  of 
the  second  impulse  the  acoustic  wave  has  set  up  an  inhomogeneous  refractive  in¬ 
dex  distribution.  If  the  distribution  is  sufficiently  converging,  the  light 
beam  from  the  second  impulse  will  be  trappc  1  at  a  small  radius. 

The  value  cf  the  model  lies  in  the  fact  that  the  sound  and  light 
equations  are  decoupled.  The  light  acts  only  once,  instantaneously,  on  the 
material.  The  resulting  sound  wave  reacts  only  once,  instantaneously,  on  the 
light  beam.  Although  this  model  does  not  represent  the  actual  phenomenon  very 
realistically,  it  is  relatively  easy  to  analyze  and  it  does  give  ^ood  value: 
for  the  trapping  threshold. 

The  model  can  be  improved  slightly  by  computing  the  sound  wave 
response  to  a  continuous  driving  force  applied  by  the  light  for  the  duration  of 
the  laser  pulse.  In  that  case,  the  sound  wave,  as  it  develops,  should  focus 
the  beam.  To  treat  the  focusing  problem,  however,  it  would  be  necessary  to 
calculate  the  radial  sound  wave  in  many  cross  sections  of  the  beam.  An  alter¬ 
native  is  to  consider  the  case  where  the  light  pulse  is  so  weak  that  the  acoustic 
waves  do  not  focus  the  beam  significantly.  It  is  then  still  possible  to  compare 
the  acoustic  response  to  a  continuous  driving  force  with  the  impulse  response. 

When  we  made  the  calculation,  we  found  that  the  response  was  basic¬ 
ally  the  same.  There  was  a  phase  delay  dua  to  the  later  arrival  of  the  majority 
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of  the  light  energy,  and  the  amplitude  of  the  response  was  reduced  slightly  as 
the  pulse  half-time  Increased.  Nevertheless,  the  same  Initial  focusing  dis¬ 
tribution  followed  by  an  outward  propagating  ripple  and  slow  settling  of  the 
wake  was  observed.  These  results  Increase  our  confidence  In  the  validity  of 
the  two-impulse  model  of  Section  III. 

The  Sound  Wave  Equation: 

Let  a  be  the  acoustic  compression,  and  r  be  the  radius  In  cylin¬ 
drical  coordinates.  The  radial  sound  wave  equation  becomes 

S2a/3r2  +  r  *50/5 r  =  d2c/d t2  -  (A  +  Bt)(l  -  r2)  exp(-r2) 

where  the  driving  term  has  the  Laplaclan  of  a  Gaussian  radial  distribution  (as 
before)  and  a  step-plus-ramp  temporal  distribution.  (The  general  pulse  tempo¬ 
ral  distribution  may  be  approximated  to  any  desired  degree  of  accuracy  by  a 
continuous  series  of  steps  and  ramps.)  The  equation  Is  Hankel  transformed 
by  multiplying  It  by  rJQ(r)  and  Integrating  from  0  to  ».  In  general,  the  Han¬ 
kel  transform  of  the  Laguerre  polynomial  Is 

J  exp(-r2)L  (r2)  rJ  (or)dr  =  (l/2n!)  (p/ 2) 2n  exp(-p2/4) 

«  no 

o 

2  2 

Our  transform  has  L^(r  )  =  1  -  r  ,  so  the  transformed  equation  becomes 
dV/dt2  +  pV  =  (A  +  Bt)  p2  exp(-p2/4)/8 

where 

o(p,")  =  J  a(r,t)r  JQ(pr)dr 
o 

Is  the  Hankel  transform  of  the  compression  and 

^*(rd2c/dr2  +  5o/dr)  J  (pr)dr  =  -p % 

i  o 

o 

Is  a  known  trar s formation.  (See  Ref.  7). 

The  transformed  equation  Is  just  the  one-dlmenslonal  wave  equation 
with  a  step-plus-ramp  driving  ‘■erm.  This  forced  equation  Is  easily  solved  by 
Laplace  transforms.  (It  could  even  be  solved  for  a  unit  quartlc  pulse  forcing 
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term.  That  solution,  however,  would  not  be  particularly  useful  later  in  the 
general  problem  of  acoustic  focusing.)  The  solution  for  the  Hankel  transform 
of  the  compression  is 

0(p,t)  =  o(p,0)  cospt  +  [dcr(p,0)/5t]  p  *  sinpt 

••1  2 

-CA(cospt-l)  +  Bp  (sinpt-ot)]  exp(-D  /^)/8 
given  the  initial  "position"  O(p,0)  and  "velocity"  C^c(p,0)/ot]  of  the  compression. 

In  order  to  evaluate  the  acoustic  response  to  a  driving  function 
of  arbitrary  temporal  distribution  by  a  continuous  step-plus-ramp  approximation, 
it  is  necessary  to  have  the  "velocity"  as  a  function  of  time  also.  Then  the 
"position"  and  "velocity"  at  the  end  of  the  timestep  become  the  new  initial 
"position"  and  "velocity"  for  the  next  timestep.  The  velocity  is  obtained  by 
simple  differentiation: 

So(p,t)/St  =  -c(p,0)  psinpt  +  [5c(p,0)/5t]  cospt 

2 

-[-Apsinpt  +  B(cospt-l)]  exp(-p  / 4 ) / 9 
Integral  Representation  of  the  Sound  Wave  Solution: 

The  acoustic  compression  is  now  obtained  by  the  inverse  Hankel 

trans  form. 

00 

r  — 

cr(r,t)  =  a(p,t)  pJ  (pr)dp 

*■'  O 

o 

oo 

=  J  ;a(p,0)  pcospt  +  [r*c(p,0)/dt]  sinpt 

o 

2 

-cAp (cospt- 1)  +  B(sinpt-pt)]  exp(-p  /4)/8j  JQ(pr)dp 

We  note  that  the  term  c(p,0)  is  multiplied  by  p  in  the  integral.  Thus,  if  Op 
is  always  evaluated,  rather  than  o  alone,  there  will  be  no  division  by  zero 
at  p=0. 

Before  the  start  of  the  optical  pulse  at  t=0,  the  medium  is  at 
rest.  Thus,  c(r,0)  =  o(p,0)  =  0.  In  th^case  of  an  extremely  short  optical 
pulse  there  may  be  an  initial  velocity  at  t=0,  as  in  Section  III.  The  initial 
velocity  due  to  an  impulse  of  Gaussian  spatially  distributed  light  is 
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da(r,0)/3t  =  (1-r2)  exp(-r2) 
with  the  Henkel  transform 

5o(p,0)/dt  =  p?  exp(-p2/4)/8 

as  before. 

Method  of  Computation: 

The  integral  representation  of  the  sound  wave  solution  above  was 
checked  by  computing  and  displaying  the  wave.  This  computer  check  was  impor¬ 
tant  because  the  integral  representation  is  a  starting  point  for  further  reduc¬ 
tions  to  the  computing  formulas  used  in  the  beam  tracing  program. 

The  only  input  variable  for  the  computation  is  the  pulse  half-time, 
measured  in  units  of  the  material  time  constant.  Time  runs  from  0  through  10 
material  time  constants  in  steps  of  0.1  time  constant.  The  unit  quartJc  pulse 
starts  from  zero  at  T=0  and  is  evaluated  after  each  timestep,  for  finite  pulse 
half-times.  If  the  pulse  half-time  is  zero,  the  program  computes  and  stores  the 
initial  velocity  da(p,0)/dt  as  above. 

For  every  pulse  halftime  entered,  the  program  computes  a  movie  of 
101  frames.  In  each  frame  the  wave  amplitude  is  plotted  at  21  points  along  the 
radius  scale,  from  0  to  5  beam  radii.  After  computation  the  maximum  amplitude 
range  was  found  to  be  -0.05  to  +0.34  for  the  impulse-excited  wave,  and  less  for 
other  excitations.  The  frame  consists  of  23  words  of  24  bits  each.  The  first 
word  contains  the  number  of  points,  not  counting  the  0th  point;  it  is  20  in  this 
case.  The  next  word  contains  the  number  of  repetitions  for  the  frame,  equal 
to  the  ratio  of  the  timestep  for  this  frame  to  the  shortest  timestep  in  the 
movie.  Since  all  the  timesteps  are  equal  in  this  program,  the  number  of  repe¬ 
titions  is  1  for  each  frame.  The  remaining  21  words  contain  the  points  packed 
for  oscilloscope  display.  The  leading  4  bits  contain  a  code  number  7,  causing 
the  oscilloscope  to  draw  intensified  vectors  from  the  last  point  displayed  to 
the  new  point.  The  next  10  bits  specify  the  scaled  value  of  the  radius,  and 
the  10  xeast  significant  bits  are  the  scaled  value  of  the  amplitude.  The  words 
are  written  in  BCD  coding,  16  octal  fields  8  characters  wide,  without  spaces. 

A  graph  of  every  4th  frame  was  also  printed  as  a  quick  check  on  the  computation. 
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The  Bessel  fu  ctlon  was  computed  according  to  a  polynomial  approxi¬ 
mation  (Ref.  4).  For  argument  Z  less  than  or  equal  to  3, 


A  =  (Z/3)2 

J  (Z)  =  1  +  A(a  +  A(a1  +  A(a  +  A(e  +  A(a  +  Aa  ))))). 
o  o  1  2  3  4  5 


For  larger  arguments, 

A  =  3/Z 


J  (Z)  =  (b  +  A(b  +  A(b.  + 

O  0  12 

•  cos(Z  +  c  +  A(c. 

o  I 


A(b3  +  A(b4  +  A(b5  +  Ab6)))))). 

+  A(c2  +  A(c3  4-  A(c^  +  A(c,_  +  Ac^))))))Z 


The  coefficients  are: 


a  =  -2.2499997 
o 

a3  =  +1.2656208 

a2  =  -0.3163866 

a3  =  +0.0444479 

a.  -0.0039444 
4 

a5  =  +0.0002100 


b  .=  +0.79788456 
o 

b  =  -0.00000077 
b2  =  -0.00552740 
b3  =  -0.000C9512 

b.  =  +0.00137273 

4 

b5  =  -0.00072805 

b,  =  +0.00014476 
6 


c  =  -0.78539816 
o 

c  3  =  -0.04166397 

c2  =  -0.00003954 

c  3  =  +0.00262573 

c  =  -0.00054125 
4 

c5  =  -0.00029333 

c .  =  +0.00013558 
6 


These  formulas  are  supposed  to  give  correct  results  to  within  ±5  x  10  absolute 


error. 


At  first  we  attempted  to  evaluate  the  integral  representation  with 
a  21-sample  Simpson  approximation.  However,  the  solution  did  not  vanish  suf¬ 
ficiently  as  T  approached  10.  After  the  ripple  passed  out  off  the  movie  to  the 
right  and  the  wake  settled,  a  low-spatial-frequeucy  ripple  seemed  to  move  back 
toward  the  left.  This  is  explained  as  insufficient  cancellation  by  the  sinpt 
and  cospt  terms,  which  have  wavelengths  of  0.628  at  T=10.  For  a  21-sample 
Simpson  approximation  over  the  range  0  to  5  the  sampling  interval  is  0.25,  or 
2.5  samples  per  wavelength.  Therefore,  we  changed  to  a  41-sample  Simpson 
approximation  over  the  same  range,  taking  5  samples  per  wavelength,  to  achieve 
satisfactory  results. 
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When  we  speak  of  a  41-sample  Simpson  approximation  we  are  including 

the  sample  at  zero,  which  does  not  contribute  to  the  i  egral.  Thus,  in  the 

entire  movie  40  samples  of  the  integral  were  ..aken  to  evaluate  each  of  21 
points  in  each  of  101  frames,  40  x  21  x  101  =  84,840  samples  in  all.  Each 
sample  requires  a  value  of  the  Bessel  function.  To  save  running  time,  we  set 
up  a  table  of  the  41  x  21  =  861  values  of  the  Bessel  function  actually  used. 

Let  i  be  the  index  for  r  values,  and  let  j  be  the  index  for  p  values.  Then 

r.  =  0. 25i,  i  =  0,1, . . . , 20 

0.  =  0.125J,  j  =  0,1,. ..,40 

Jij  =  Jo(ri°j)  =  VO-031251  j) 

Clearly,  =  J^,  and  =  80  the  table  contains  many  duplicate 

entries.  However,  we  could  not  think  of  a  simple  indexing  scheme  to  reduce 
the  number  of  entries  by  eliminating  redundancy,  so  we  stored  the  entire  table 
as  it  stands  in  computer  memory  at  the  start  of  the  program.  Various  compli¬ 
cated  indexing  schemes  are  possible,  of  course,  but  none  of  the  ones  we  thought 
of  could  be  computed  qulrkly. 

Three  additional  working  arrays  are  required  to  compute  the  integral. 
Each  array  has  41  elements.  It  is  necessary  to  save  the  values  of  pcr(p,  t  = 
timestep)  and  ?ij(p,t=timestep) />t  after  each  frame,  in  order  to  have  prf(o.O) 
and  da( p,0) At  for  the  next  frame.  Also,  the  kernel  of  the  integral  (apart 
from  the  factor  Jq  (pr))  does  not  involve  r,  so  41  samples  of  it  may  be  saved 
and  used  repeatedly  in  each  frame.  It  is  useful  to  include  the  Simpson  weight¬ 
ing  coefficients  as  part  of  the  kernel. 

Finally,  the  21  points  in  each  frame  were  stored  so  that  they  could  be 
written  as  a  single  record  for  each  frame. 

The  program  is  carried  out  i.i  steps  as  follows.  Here  the  equal  sign 
means  place  the  value  of  the  expression  or.  the  right  into  the  cell  named  on  the 
left. 

1.  Create  the  Bessel  function  table. 

Jjj  =  (ijdpdr),  i  =  0,1,. ..,20;  j  =  0,1,. ..,40 

where  dp  =  1/8  and  dr  =  1/4. 
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2.  Read  p,  the  half-time  of  the  pulse.  If  there  are  no  more 
values  of  p  to  read,  stop. 

3.  Clear  working  arrays. 

t  =  0;  h.  =  kj  =  0,  j  =  0,1,.... .40 
where  t  is  the  time,  h^  is  a(p^,0) 
and  kj  is  the  jth  sample  of  the  kernel. 

4.  Set  the  initial  'velocity". 

d.  =0  if  the  p  is  non-zero 
J  2  2 

d.  =  p  exp(-p ./4)/8,  j  =  0,1,. ..,40  is  p  is  zero 


where  d.  is  dc(p. , 0)/5t 
J 


5.  Go  to  9. 

6.  Compute  the  step  and  ramp  coefficients, 
a  =  q(t-dt;p),  b  =  Cq(t ;p)-a]/dt 

where  q  is  the  unit  quartic  pulse  of  Section  IV,  and  at 
is  the  timestep,  1/10. 

7.  Compute  the  kernel.  Also,  compute  the  values  of  "position" 
and  "velocity"  required  for  the  next  timestep. 

Do  the  next  6  replacements  successively  as  j  takes  on  the 
values  1, 2, ... ,40. 

e  =  exp(-Pj/4)/8 

f  =  apj (cosp^dt-i)  +  b(sinPjdt-Pjdt) 

g  =  h.cosp.dt  +  d.sinp.dt  -  ef 
J  j  J  J 

kj  =  (1/3)  dp  2(jmod2+l)g 

d.  =  -h.sinp.dt  +  d.cosp.dt  -  e[-an ,sinp  ,dt  +  b(cosp,dt  -  1)] 
J  j  j  J  j  J  j 

hj  '6 

Finally,  adjust  the  end  sample  of  the  kernel. 
k40  *  k40/2 


8.  Compute  the  sound  wave  solution  at  21  radial  points. 
Simultaneous!)  print  them  every  4th  timestep. 
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Si=Sj=o  Jljkj’  i.-  0,1,. ..,40 


9.  Pack  the  points  as  described  before  Into  21  words  and  write 
them  on  magnetic  tape. 

10.  Take  the  next  tlmestep.  If  this  Is  the  last  frame,  go  to  2. 
Otherwise, 
t  =  t  +  dt 
and  go  to  6. 


This  87-statement  SDS  Fortran  IV  program  required  1  minute  54 
seconds  for  compilation  and  4  minutes  12  seconds  to  load  and  run  four  cases 
with  p  equal  to  0,  0.5,  1,  and  1.5.  Later  we  used  a  second  program  to  place 
corresponding  frames  of  the  four  movies  In  four  quadrants  of  a  single  frame 
of  a  new  movie.  The  latter  movie  permitted  us  to  make  phase  and  amplitude 
comparisons  for  the  four  cases. 


The  results  are  plotted  and  reproduced  In  Figures  19  through  24, 
for  t  =  0  to  t  =  7. 1.  The  upper  left  quadrant  of  each  frame  shows  the  response 
for  p  -  0,  the  upper  left  for  p  =  0.5,  the  lo^er  left  for  p  =  1,  and  the  lower 
right  for  p  =  1.5. 

6.2  THE  LOCAL  LENS  ESTABLISHED  BY  A  SINGLE  LASER  IMPULSE 

A  laser  pulse  whose  half-time  v;as  negligible  compared  with  the 
material  time  constant  would  pass  through  a  transparent  dielectric  on  a  hyper¬ 
bolic  path  without  being  focused  appreciably  by  electrostrlctlon,  because  of 
the  time  required  for  acoustic  compression.  Nevertheless,  the  medium  would 
be  set  Into  motion,  and  the  focusing  power  would  develop  with  time.  An  analysis 
of  this  model  shows  how  a  single  Impulse  can  set  up  e  localized  acoustic  ler.s 
at  the  Initial  focus  of  the  laser  beam.  Because  the  acoustic  lens  Is  localized, 
It  can  cause  a  second  focus  for  the  laser  beam,  provided  that  the  next  laser 
Impulse  comes  before  the  sound  wave  has  subsided.  This  behavior  may  be  con¬ 
trasted  with  the  behavior  In  Kerr  effect  trapping,  where  only  one  focus  would 
occur  under  such  conditions. 


The  second  term  on  the  right-hand  side  of  equation  (2)  of  Section 
IV  Is  the  focusing  term  In  the  beam  tracing  equation, 
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Figure  24.  Amplitude  Versus  Radius  For 
Four  Sound  Waves,  T  =  6.0  to 
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(a/n^)  52n(0, z,t)/5r" 

Derivatives  of  the  refractive  index  are  proportional  to  the  acoustic  compression. 

These  can  be  evaluated  analytically  at  r  =  0.  The  result  will  be  the  dependence 
of  the  focusing  term  on  the  initial  hyperbolic  beam  radius. 

I 

We  are  therefore  interested  in  the  second  derivative  of  the  acous¬ 
tic  compression  along  the  beam  axis  as  a  function  of  z  and  t.  This  quantity 
is  proportional  to  the  coefficient  of  the  focusing  term.  We  will  use  the 
integral  representation  of  the  compression  when  excited  only  by  an  initial 

impulse.  To  get  the  full  dependence  on  beam  radius,  we  will  work  the  problem  I 

through  in  unsealed,  physical  coordinates. 

The  path  taken  by  the  beam  through  a  homogeneous  medium  is 

a(z)  =  C(a  +  a  z)2  +  (\z/2nn  a  )2]^2 
o  i  o  o 

as  in  equation  (3),  Section  IV.  The  laser  beam  intensity  is 

I(r,z,t)  =■•  (W/ita2)  exp(-r2/a2)  [6(t)  +  fl(t-5p/8)]/2 
where  the  spacing  5p/8  between  delta  functions  corresponds  to  placing  an  impulse 
at  the  centroid  of  each  half  of  the  unit  quartic  pulse. 

The  force  exerted  on  the  medium  by  elt ctrostriction,  assuming  the 
Clausius  -  Mosotti  relation  applies,  is 

f  =  (n2  +  2)  (n2  -  1)  V  I/3n  c 
o  o  o 

The  velocity  field  after  the  first  impulse  is 

o+ 

r*  • 

j  f  dt  =  g  u  (r,z,0+) 

where  g  is  the  average  density  of  the  medium,  by  conservation  of  linear  momen¬ 
tum,  and  the  compression  is  the  negative  divergence  of  the  displacement, 

— *  — ♦ 

C  =  -  V  •  u 

Thus,  in  cylindrical  coordinates,  taking  the  time  derivative  of  the  above 
equation, 
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dc(r,z,0+)/dt  =  -  V  •  u  (r,z,0+)  q+ 

=  -(n  +  2)(n  -1)  V2  j  x  dt/3  n  eg 

0  °  .1  o 

=  -[W(n2  +  2)(n2  -  l)/6n  cgna2] 
o  o  o 

{r  *  d[rd  exp(-r2/a")/eJr]/5rj 
-  (2W/3jt  nocga4)  (n2  +  2)  (n2  -  1)  (1  -  r2/a2)exp(-r2/a2) 

The  Hankel  transform  Is  easily  evaluated  as 

da(p,z,0+)/dt  =  (2W/3«n  ega2)  (n2  +  2)  (n2  -  1)  • 

o  o  c 

pco 

j  (r/a)  Jq[ (ap) (r/ a) ]  (1  -  r2/a2)  exp(-r2/a2)  d(r/a) 
o 

=  (W/12  it  nQCg)  (n2  +  2)  (n2  -  1)  p2  exp(-a2p2/4) 

The  material  was  In  a  state  of  rest  before  the  initial  impulse,  so  0(p,z,C+)  0. 

Also,  there  Is  no  forcing  term  after  the  impulse,  as  the  intensity  drops  to  zero 
again,  so  A  =  3  =  0.  Thus,  by  substitution  in  the  integral  representation  of 
the  solution  we  obtain 

a(r,z,t)  =  j  [dr7(p,*,0+)/dt]  sin(pvt)  J  (pr)  dp 
o 

CD 

"  2  2  2 

=  c  j  p  exp  ("S'  p  /4)  sln(pvt)  J  (pr)  dp 
o  ° 

where  we  have  inserted  the  velocity  of  sound  v  and 

C  =  (W/I2nn  cgv)  (n2  +  2)  (n2  -  1) 
o  o  o 

2  2 

We  may  now  obtain  *  a(0,x,t)/dr  by  differentiating  under  the  integral  sign. 

We  have 

92J  (pr)/?^r2  =  [j„  (pr)  -  J  (pr)]  p2/2 

O  L  O 

and 

J2(0)  =  0,  JQ(0)  =  1. 
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This  leads  to 

09 

d2fT(,?-,z,t)/dr2  -  -(C/2)  J  exp(-a2pV4)  sin(pvt)dp 

o 

The  remaining  integral  is  just  a  Fourier  transform.  A  useful  property  of 
Fourier  transforms  permits  us  tc  simplify  the  integral  further.  If 

F(y)  «■  J  f (x)  sin(xy)  dx  for  y  >  0 
o 

is  the  Fourier  transform  of  f(x),  then  the  Fourier  transform  of  x2m  f(x)  is 
(-l)m  d2m  F(y)/dy2m.  Thus  we  need  or  ly  evaluate 


K  = 


09  03 

2  2  f  2  2 

exp(~a  p  / 4)  sin(pvt)  dp  =  Im  ]  exp(-a  o  /4  +  iovt)dp 

L  «J  . 

O 


where  ImL  ]  denotes  the  imaginary  part.  By  completing  the  square  in  the 
exponential  and  shifting  the  variable  we  can  reduce  the  remaining  part  of 
the  integral  to  a  finite  interval. 


K  =  Im  J  exp[-  (ap/2  -  ivt/a)2J  dp  exp(-v2t2/a2) 1 


Let  P  =  ap/2  -  ivt/a;  then 


K  Im  -^(2/a)  exp(-v2t2/a 2)  ,  _ivt/a  exP("^)d£j 


-ivt/a 


=  Im  '^(2/a'  exp(-v2t2/a2)^J  exp(-"2)d!r  -  exp(-F2)dp 


Now  let  “  -§/i  and  take  only  the  imaginary  part 

->vt/a 


2  2  2 

K  =  (2/'a)  exp(-v  t  /a  )  j  exp(rf)  df] 

o 


With  the  above  result,  the  focusing  term  becomes 
5  0(O,z,t)/dr2  -  -  (C,/2)  d  *K/d(vt)^ 
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It  is  now  convenient  to  change  .o  scaled  coordinates.  We  l«it  T  =  vt/a  so 

*r 

2  2  4  p  2  r*  2  "  a 

9  a(0,z,t)/?*r  =  -Ca  d  ^exp(-T  )  exp(r  )dri  ./dT 

o 


=  -Ca"5  20t 

U 


8t3  +  (12  -  48t2  +  16T4)exp(-T2)j  expCn^dr) 

o 


For  very  small  values  of  T  we  may  replace  the  exponentials  by  1  and  drop  all 
terms  of  order  T2  or  less. 

Then, 


92c(0,z,t  «  a/v)/5r2 


-  5  -6 

■8Ca  t  =  ~8Cvt  a 


As  the  terms  in  higher  powers  of  t  become  dominant,  the  coefficient  of  the 
focusing  term  becomes  even  more  localized  near  the  original  focus. 

2  rT  2 

The  Dawson  integral  D(t)  =  exp(-T  )  J  exp(Ti  )dt)  is  tabulated  (Ref.  4, 
p.  319).  For  reference,  other  derivatives  are  a§  follows: 

dD/dT  •—  1  —  2tD 
d2D/dT2=  -2t  +  (4t2  -  2)D 
d3D/dT3=  4t^  -  4  +  (12t  -  8t3)D 
d4D/dT4-  20t  -  8t3  +  (12  -  48t2  +  16t4)D 

An  asymptotic  expansion  for  D  is 

2 


Lim  D(t  -  *)  =  -i  exp(-T  )  /jT/ 2  +  (1/2t)  [l  +  1/2t  +  1 

+  1 


3/(2T2)2 


3  •  5/ (2t2)3  h...] 


Each  of  the  derivatives  of  D  becomes  0  as  t  approaches  infinity  because  the 
polynomial  terms  are  cancelled  bv  subtraction.  Thus,  the  formulas  are  mathe¬ 
matically  correct,  but  inconvenient  for  computation  for  large  values  of  t. 

The  coefficient  of  the  focusing  term  was  computed  and  plotted 
versus  z  in  Figure  25.  The  function  of  a  was  given  by  equation  (3)  of  Section 
IV,  plotted  in  the  first  frame  of  Figure  14.  In  Figure  25,  reading  left  to 
right  and  top  to  bottom,  we  see  the  local  development  of  the  focusing  coeffi¬ 
cient  as  t  runs  from  0  to  2.  In  the  two-impulse  model  of  trapping,  the  second 
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Figure  25.  Coefficient  of  the  Focusing  Term  Plotted  Versus  Axial  Length 
At  Various  Fractions  ;  f  the  Material  Time  Constant,  After  an 
Initial  Impulse 
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impulse  could  arrive  anytime  before  t  =  12  and  still  pass  through  a  converging 
acoustic  lens.  If  the  pulse  half-time  were  much  longer,  it  would  first  be 
focused,  then  defocused- 

In  the  plot  the  coefficient  at  the  initial  focus  appears  to 
develop  with  a  phase  just  twice  as  fast  as  the  coefficient  at  the  left-hand 
edge,  and  the  amplitude  it  32  times  greater.  This  is  because  of  the  2  to  1 
ratio  of  values  of  a. 

The  Dawson  in'egral  v;as  computed  by  running  Simpson  integration, 
using  a  step  size  of  0.0005.  (A  step  size  of  0.05  was  too  large  for  stable 
computation  beyond  about  T  =  4. ) 

Let 

0.001J 

d ,  ~  J  expCn2)dr| 
j  o 

Than 

d  =  0 
o 

t  =  0.001J  d  =  d  +  (0.0005/3)  |exp[ (t-0.001)2]  +  4  rxpZ (t-0. 00  '5) 2]  +  ex 

defines  the  running  Simpson  integral  used. 

6.3  STEP-BK-STEP  SOLUTION  OF  THE  COUPLED  EQUATIONS 

The  solution  of  the  coupled  sound  and  beam  radius  equations  has  al¬ 
ready  been  discussed  in  Section  V.  The  actual  formulas  used  by  the  computer 
program  will  be  presented  here. 

The  step-plus -ramp  driven  sound  wave  solution  required  becomes 

7  2  a® 

9  <?( 0,z,t)/dr  =  —  ( 1/2)  j  p2  {fl*(p,z,0)pcospt  +  [dc(p, z,0)/dt]  sinpt 

o 

2  2  2  'I 

-(a  /8)  exp(-a  o  /4)  [Aap(cos  aot-1)  +  B(sin  apt-apt) j [dp 

o(p,z,t)p  =  7(p,z,0,p  cospt  +  [3a(p,z,0)/9t]  sin  pt 

2  2  2 

-(a  / 8)  exp(-a  p  /4)[Aap(c  os  spt-1)  +  B(sin  aot-apt)] 
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&j(p,z,t)/5t  =  -o(p, z,0)p  sin  pt  +  [?^r(p,z,  0)/?)t  J  cos  pt 

■(a^/8)  exp(-a^p^/4) r"Aap  sin  apt  +  B(cos  apt  -1)] 

These  may  be  evaluated  by  Simpson  integration  as  Jescribed  before  in  Subseccion 
6.1.  The  equations  now  include  a  =  a(z,t),  which  is  obi  aired  by  beam  tracing 
before  each  timestep.  The  coefficient  of  electrostriction  can  be  entered  sep¬ 
arately,  because  the  sound  wave  equation  is  linear. 

For  the  Runge-Kutt  solution  of  the  beam  tracing  equation  it  was 

necessary  to  interpolate  between  points  where  the  focusing  term  was  known.  Ue 

used  the  Lagrange  cubic  fitting  formula.  If  interpolation  is  required  between 

X  and  X  . ,  where  Y  ,,  Y  ,  Y  ,  Y  are  known  ordinates  for  equally  spaced 
n  n+1  n- 1  n  n+ij  n-r^ 

abscissas,  and  £  =  (X-X  Y/h  is  the  fraction  of  the  distance  from  the  last  abscissa 
to  the  point  of  interpolation  divided  by  the  step  size, 

then, 

Y(Xn  +  ?h)  =  -P(ff-l)(f-2)  Y^/6  +  (~2“1) (?”2)  Yr/2 
-f<*+l)(?-2)  Yn+1/2  +  ?(^2-l)  Yn+2/6 

See  Ref.  4,  pp.  878f. 
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